# for a cython tutorial, see: 
# http://conference.scipy.org/proceedings/SciPy2009/paper_1/full_text.pdf

cdef extern from "math.h":
    double exp(double)
    double cos(double)
    double sin(double)
    double abs(double)
    double sqrt(double)
    double pow(double base, double exponent)

import numpy as np
cimport numpy as np

cdef double pi = np.pi

# sign function
# =============

cdef int sign(double x):
    if x > 0.:
        return 1
    elif x < 0.:
        return -1
    else:
        return 0
  


# cutoff function and derivative
# ==============================

cdef double f_cut(double x, double cutoff, double d) except *:
    if x <= cutoff - d:
        return 1.
    elif cutoff - d < x <= cutoff:
        return (cos(pi * (x-cutoff+d)/d) + 1.)/2.
    else:
        return 0.

cdef double dx_f_cut(double x, double cutoff, double d) except *:
    if cutoff - d < x <= cutoff:
        return -pi/(2.*d) * sin(pi * (x-cutoff+d)/d)
    else:
        return 0.

# squared exponential (SE) covariance function and derivatives
# ============================================================

cdef double k_SE(double x1, double x2, double delta,
                 double theta) except *:
    return delta**2*exp(-(x1-x2)**2/(theta)**2)

cdef double dx1_k_SE(double x1, double x2, double delta,
                     double theta) except *:
    return (x2-x1)/theta**2 * k_SE(x1, x2, delta, theta)

cdef double dx1_dx2_k_SE(double x1, double x2, double delta,
                         double theta) except *:
    return (theta**2 - (x1-x2)**2/theta**4) * k_SE(x1, x2, delta, theta)

# damped squared exponential (SE_damped) covariance function and derivatives
# ==========================================================================

cdef double k_SE_damped(double x1, double x2, double delta,
                        double theta, double beta) except *:
    return exp(-beta*(x1+x2))*delta**2*exp(-(x1-x2)**2/(theta)**2)

cdef double dx1_k_SE_damped(double x1, double x2, double delta,
                            double theta, double beta) except *:
    return ((x2-x1)/theta**2 - beta) * k_SE_damped(x1, x2, delta, theta, beta)

cdef double dx1_dx2_k_SE_damped(double x1, double x2, double delta,
                                double theta, double beta) except *:
    return (1./theta**2 - ((x2-x1)/theta**2 - beta)**2)*k_SE_damped(x1, x2, 
                                                            delta, theta, beta)

# laplacian covariance function and derivatives
# =============================================

cdef double k_Lap(double x1, double x2, double delta,
                  double theta) except *:
    return delta**2 * exp(-sqrt((x1-x2)**2)/theta)

cdef double dx1_k_Lap(double x1, double x2, double delta,
                      double theta) except *:
    return sign(x2-x1)/theta * k_Lap(x1, x2, delta, theta)

cdef double dx1_dx2_k_Lap(double x1, double x2, double delta,
                          double theta) except *:
    return -1./theta**2 * k_Lap(x2, x2, delta, theta)

# rational quadratic (RQ) covariance function and derivatives
# ===========================================================

cdef double k_RQ(double x1, double x2, double delta,
                 double theta, double beta) except *:
    return delta**2 * pow((1. + (x1-x2)**2/(2*beta*theta**2)), -beta)

cdef double dx1_k_RQ(double x1, double x2, double delta,
                     double theta, double beta) except *:
    return (x2-x1)/theta**2 * pow((1. + (x1-x2)**2/(2*beta*theta**2)), -(beta+1.))

cdef double dx1_dx2_k_RQ(double x1, double x2, double delta,
                         double theta, double beta) except *:
    return (1./theta**2 * pow((1. + (x1-x2)**2/(2*beta*theta**2)), -(beta+1.))
            - (beta+1.)/beta * (x1-x2)**2/theta**4
            * pow((1. + (x1-x2)**2/(2*beta*theta**2)), -(beta+2.)))


# definition of kernel objects for use in covariance matrix build functions
# =========================================================================

# base class for kernel functions
cdef class Kernel:
    cpdef double evaluate(self, 
                          double x1, 
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return 0

cdef class K_SE(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return k_SE(x1, x2, delta, theta)

cdef class dx1_K_SE(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return dx1_k_SE(x1, x2, delta, theta)

cdef class dx1_dx2_K_SE(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return dx1_dx2_k_SE(x1, x2, delta, theta)

cdef class K_SE_cut(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return (k_SE(x1, x2, delta, theta)
                * f_cut(x1, cutoff, d)
                * f_cut(x2, cutoff, d))

cdef class dx1_K_SE_cut(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return (dx1_k_SE(x1, x2, delta, theta)
                * f_cut(x1, cutoff, d) * f_cut(x2, cutoff, d)
                + k_SE(x1, x2, delta, theta)
                * dx_f_cut(x1, cutoff, d) * f_cut(x2, cutoff, d))

cdef class dx1_dx2_K_SE_cut(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return (dx1_dx2_k_SE(x1, x2, delta, theta)
                * f_cut(x1, cutoff, d) * f_cut(x2, cutoff, d)
                + dx1_k_SE(x1, x2, delta, theta)
                * f_cut(x1, cutoff, d) * dx_f_cut(x2, cutoff, d)
                - dx1_k_SE(x1, x2, delta, theta)
                * dx_f_cut(x1, cutoff, d) * f_cut(x2, cutoff, d)
                + k_SE(x1, x2, delta, theta)
                * dx_f_cut(x1, cutoff, d) * dx_f_cut(x2, cutoff, d))

cdef class K_SE_damped(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return k_SE_damped(x1, x2, delta, theta, beta)

cdef class dx1_K_SE_damped(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return dx1_k_SE_damped(x1, x2, delta, theta, beta)

cdef class dx1_dx2_K_SE_damped(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return dx1_dx2_k_SE_damped(x1, x2, delta, theta, beta)

cdef class K_SE_damped_cut(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return (k_SE_damped(x1, x2, delta, theta, beta)
                * f_cut(x1, cutoff, d)
                * f_cut(x2, cutoff, d))

cdef class dx1_K_SE_damped_cut(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return (dx1_k_SE_damped(x1, x2, delta, theta, beta)
                * f_cut(x1, cutoff, d) * f_cut(x2, cutoff, d)
                + k_SE_damped(x1, x2, delta, theta, beta)
                * dx_f_cut(x1, cutoff, d) * f_cut(x2, cutoff, d))

cdef class dx1_dx2_K_SE_damped_cut(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return (dx1_dx2_k_SE_damped(x1, x2, delta, theta, beta)
                * f_cut(x1, cutoff, d) * f_cut(x2, cutoff, d)
                + dx1_k_SE_damped(x1, x2, delta, theta, beta)
                * f_cut(x1, cutoff, d) * dx_f_cut(x2, cutoff, d)
                - dx1_k_SE_damped(x1, x2, delta, theta, beta)
                * dx_f_cut(x1, cutoff, d) * f_cut(x2, cutoff, d)
                + k_SE_damped(x1, x2, delta, theta, beta)
                * dx_f_cut(x1, cutoff, d) * dx_f_cut(x2, cutoff, d))

cdef class K_Lap(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return k_Lap(x1, x2, delta, theta)

cdef class dx1_K_Lap(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return dx1_k_Lap(x1, x2, delta, theta)

cdef class dx1_dx2_K_Lap(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return dx1_dx2_k_Lap(x1, x2, delta, theta)

cdef class K_Lap_cut(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return (k_Lap(x1, x2, delta, theta)
                * f_cut(x1, cutoff, d)
                * f_cut(x2, cutoff, d))

cdef class dx1_K_Lap_cut(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return (dx1_k_Lap(x1, x2, delta, theta)
                * f_cut(x1, cutoff, d) * f_cut(x2, cutoff, d)
                + k_Lap(x1, x2, delta, theta)
                * dx_f_cut(x1, cutoff, d) * f_cut(x2, cutoff, d))

cdef class dx1_dx2_K_Lap_cut(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return (dx1_dx2_k_Lap(x1, x2, delta, theta)
                * f_cut(x1, cutoff, d) * f_cut(x2, cutoff, d)
                + dx1_k_Lap(x1, x2, delta, theta)
                * f_cut(x1, cutoff, d) * dx_f_cut(x2, cutoff, d)
                - dx1_k_Lap(x1, x2, delta, theta)
                * dx_f_cut(x1, cutoff, d) * f_cut(x2, cutoff, d)
                + k_Lap(x1, x2, delta, theta)
                * dx_f_cut(x1, cutoff, d) * dx_f_cut(x2, cutoff, d))

cdef class K_RQ(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return k_RQ(x1, x2, delta, theta, beta)

cdef class dx1_K_RQ(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return dx1_k_RQ(x1, x2, delta, theta, beta)

cdef class dx1_dx2_K_RQ(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return dx1_dx2_k_RQ(x1, x2, delta, theta, beta)

cdef class K_RQ_cut(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return (k_RQ(x1, x2, delta, theta, beta)
                * f_cut(x1, cutoff, d)
                * f_cut(x2, cutoff, d))

cdef class dx1_K_RQ_cut(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return (dx1_k_RQ(x1, x2, delta, theta, beta)
                * f_cut(x1, cutoff, d) * f_cut(x2, cutoff, d)
                + k_RQ(x1, x2, delta, theta, beta)
                * dx_f_cut(x1, cutoff, d) * f_cut(x2, cutoff, d))

cdef class dx1_dx2_K_RQ_cut(Kernel):
    cpdef double evaluate(self,
                          double x1,
                          double x2,
                          double delta,
                          double theta,
                          double beta,
                          double cutoff,
                          double d) except *:
        return (dx1_dx2_k_RQ(x1, x2, delta, theta, beta)
                * f_cut(x1, cutoff, d) * f_cut(x2, cutoff, d)
                + dx1_k_RQ(x1, x2, delta, theta, beta)
                * f_cut(x1, cutoff, d) * dx_f_cut(x2, cutoff, d)
                - dx1_k_RQ(x1, x2, delta, theta, beta)
                * dx_f_cut(x1, cutoff, d) * f_cut(x2, cutoff, d)
                + k_RQ(x1, x2, delta, theta, beta)
                * dx_f_cut(x1, cutoff, d) * dx_f_cut(x2, cutoff, d))

# build K_M matrix
# ================

def calc_K_M(Kernel k,
             bond_types,
             bond_dict,
             double epsilon,
             int M):
    cdef np.ndarray[np.float64_t, ndim=2] K_M = np.zeros((M,M), 
                                                         dtype=np.float64)
    cdef int M_bond
    cdef int offset = 0
    cdef np.ndarray[np.float64_t, ndim=1] X_pseudo_bond
    cdef double delta_bond
    cdef double theta_bond
    cdef double beta_bond
    cdef double cutoff_bond
    cdef double d_bond
    cdef Py_ssize_t i, j
    # iterate over every bond type
    for b in bond_types:
        M_bond = bond_dict[b]['M']
        X_pseudo_bond = bond_dict[b]['X_pseudo']
        delta_bond = bond_dict[b]['delta']
        theta_bond = bond_dict[b]['theta']
        beta_bond = bond_dict[b]['beta']
        cutoff_bond = bond_dict[b]['cutoff']
        d_bond = bond_dict[b]['d']
        # construct covariance submatrix for current bond type
        for i in range(M_bond):
            X_pseudo_bond_i = X_pseudo_bond[i]
            K_M[i+offset,i+offset] = k.evaluate(X_pseudo_bond_i,
                                                X_pseudo_bond_i,
                                                delta_bond,
                                                theta_bond,
                                                beta_bond,
                                                cutoff_bond,
                                                d_bond) + epsilon
            for j in range(i+1, M_bond):
                K_M[i+offset,j+offset] = k.evaluate(X_pseudo_bond_i,
                                                    X_pseudo_bond[j],
                                                    delta_bond,
                                                    theta_bond,
                                                    beta_bond,
                                                    cutoff_bond,
                                                    d_bond)
                K_M[j+offset,i+offset] = K_M[i+offset,j+offset]
        
        offset += M_bond
    
    return K_M

        
# build K_NM matrix
# =================

def calc_K_NM_1bF(Kernel k,
                  bond_types,
                  bond_dict,
                  int N,
                  int M,
                  B,
                  X):
    cdef np.ndarray[np.float64_t, ndim=2] K_NM = \
            np.zeros((N,M), dtype=np.float64)
    cdef np.ndarray[np.float64_t, ndim=1] X_pseudo_bond
    cdef int i, iR, L, offset
    cdef Py_ssize_t j, l
    cdef np.ndarray[np.float64_t, ndim=1] Xi
    #cdef np.ndarray[np.float64_t, ndim=1] Di
    # iterate over atoms 
    for i in range(N):
        offset = 0
        #iR = (i-i%3)/3
        Bi = B[i]
        Xi = X[i]
        #Di = D[i]
        L = Xi.shape[0]
        # iterate over atom types of pseudo inputs
        for b in bond_types:
            # unpacking hyperparameters and pseudo inputs
            M_bond = bond_dict[b]['M']
            X_pseudo_bond = bond_dict[b]['X_pseudo']
            delta_bond = bond_dict[b]['delta']
            theta_bond = bond_dict[b]['theta']
            beta_bond = bond_dict[b]['beta']
            cutoff_bond = bond_dict[b]['cutoff']
            d_bond = bond_dict[b]['d']

            # iterate over pseudo inputs of current bond type
            for j in range(M_bond):
                X_pseudo_bond_j = X_pseudo_bond[j]
                for l in range(L):
                    # skip nonequal atom-types (they have zero covariance)
                    if Bi[l] != b:
                        continue
                    # calculate covariance
                    K_NM[i,j+offset] +=         k.evaluate(Xi[l],
                                                       X_pseudo_bond_j,
                                                       delta_bond,
                                                       theta_bond,
                                                       beta_bond,
                                                       cutoff_bond,
                                                       d_bond)
            offset += M_bond
    return K_NM


def calc_K_NM_2bF(Kernel k, 
                  bond_types, 
                  bond_dict,
                  int N,
                  int M,
                  B,
                  X,
                  D):
    cdef np.ndarray[np.float64_t, ndim=2] K_NM = \
            np.zeros((N,M), dtype=np.float64)
    cdef np.ndarray[np.float64_t, ndim=1] X_pseudo_bond
    cdef int i, iR, L, offset
    cdef Py_ssize_t j, l
    cdef np.ndarray[np.float64_t, ndim=1] Xi
    cdef np.ndarray[np.float64_t, ndim=1] Di
# iterate over total force components
    for i in range(N):
        offset = 0
        iR = int((i-i%3)/3)
        Bi = B[iR]
        Xi = X[iR]
        Di = D[i]
        L = Di.shape[0]
        # iterate over bond types of pseudo inputs
        for b in bond_types:
            # unpacking hyperparameters and pseudo inputs
            M_bond = bond_dict[b]['M']
            X_pseudo_bond = bond_dict[b]['X_pseudo']
            delta_bond = bond_dict[b]['delta']
            theta_bond = bond_dict[b]['theta']
            beta_bond = bond_dict[b]['beta']
            cutoff_bond = bond_dict[b]['cutoff']
            d_bond = bond_dict[b]['d']

            # iterate over pseudo inputs of current bond type
            for j in range(M_bond):
                X_pseudo_bond_j = X_pseudo_bond[j]
                for l in range(L):
                    # skip nonequal bond-types (they have zero covariance)
                    if Bi[l] != b:
                        continue
                    # calculate covariance
                    K_NM[i,j+offset] += Di[l] * k.evaluate(Xi[l],
                                                           X_pseudo_bond_j,
                                                           delta_bond,
                                                           theta_bond,
                                                           beta_bond,
                                                           cutoff_bond,
                                                           d_bond)
            offset += M_bond
    return K_NM


# prediction
# ==========

def vpredict(Kernel k,
             bond,
             bond_types,
             bond_dict,
             int M,
             np.ndarray[np.float64_t, ndim=1] alpha,
             np.ndarray[np.float64_t, ndim=1] X_star):
    
    cdef np.ndarray[np.float64_t, ndim=1] k_star = np.zeros(M, dtype=np.float64)
    cdef int P = X_star.shape[0]
    cdef np.ndarray[np.float64_t, ndim=1] result = np.zeros(P, dtype=np.float64)
    cdef np.ndarray[np.float64_t, ndim=1] X_pseudo_bond
    cdef Py_ssize_t i, j
    
    # determine start index for the M_bond non-zero elements of k_star
    index_of_bond = bond_types.index(bond)
    offset = 0
    for i in range(index_of_bond):
        offset += bond_dict[bond_types[i]]['M']

    # unpacking hyperparameters and pseudo inputs
    M_bond = bond_dict[bond]['M']
    X_pseudo_bond = bond_dict[bond]['X_pseudo']
    delta_bond = bond_dict[bond]['delta']
    theta_bond = bond_dict[bond]['theta']
    beta_bond = bond_dict[bond]['beta']
    cutoff_bond = bond_dict[bond]['cutoff']
    d_bond = bond_dict[bond]['d']
    
    # calculate result
    for i in range(P):
        xi = X_star[i]
        for j in range(M_bond):
            k_star[j+offset] = k.evaluate(xi, X_pseudo_bond[j],
                                          delta_bond, 
                                          theta_bond, 
                                          beta_bond, 
                                          cutoff_bond, 
                                          d_bond)
        result[i] = np.dot(k_star, alpha)
    
    return result

def system_2bF_predict(Kernel k,
                       bond_types,
                       bond_dict,
                       int M,
                       np.ndarray[np.float64_t, ndim=1] alpha,
                       B_val,
                       X_val,
                       D_val):
    cdef int N_val = len(D_val)
    cdef int i, iR, L, offset
    cdef np.ndarray[np.float64_t, ndim=1] X_val_i, D_val_i, X_pseudo_bond
    cdef Py_ssize_t j, l
    cdef np.ndarray[np.float64_t, ndim=1] k_star = np.zeros(M, dtype=np.float64)
    cdef np.ndarray[np.float64_t, ndim=1] result = np.zeros(N_val, dtype=np.float64)
    cdef double temp
    for i in range(N_val):
        offset = 0
        iR = int((i-i%3)/3)
        B_val_i = B_val[iR]
        X_val_i = X_val[iR]
        D_val_i = D_val[i]
        L = D_val_i.shape[0]
        # iterate over bond types 
        for b in bond_types:
            # unpacking hyperparameters and pseudo inputs
            M_bond = bond_dict[b]['M']
            X_pseudo_bond = bond_dict[b]['X_pseudo']
            delta_bond = bond_dict[b]['delta']
            theta_bond = bond_dict[b]['theta']
            beta_bond = bond_dict[b]['beta']
            cutoff_bond = bond_dict[b]['cutoff']
            d_bond = bond_dict[b]['d']
            # iterate over pseudo inputs of current bond type
            for j in range(M_bond):
                X_pseudo_bond_j = X_pseudo_bond[j]
                temp = 0.
                for l in range(L):
                    # skip nonequal bond-types (they have zero covariance)
                    if B_val_i[l] != b:
                        continue
                    temp += D_val_i[l] * k.evaluate(X_val_i[l], X_pseudo_bond_j,
                                                    delta_bond, 
                                                    theta_bond, 
                                                    beta_bond, 
                                                    cutoff_bond, 
                                                    d_bond)
                k_star[j+offset] = temp
            
            offset += M_bond
        
        result[i] = np.dot(k_star, alpha)
    
    return result

def system_1bE_predict(Kernel k,
                       bond_types,
                       bond_dict,
                       int M,
                       np.ndarray[np.float64_t, ndim=1] alpha,
                       B_val,
                       X_val):
    cdef int i, iR, L, offset
    cdef np.ndarray[np.float64_t, ndim=1] X_val_i, X_pseudo_bond
    cdef Py_ssize_t j, l
    cdef np.ndarray[np.float64_t, ndim=1] k_star = np.zeros(M, dtype=np.float64)
    #cdef np.ndarray[np.float64_t, ndim=1] result = np.zeros(1, dtype=np.float64)
    cdef double temp
 
    offset = 0
    B_val_i = B_val
    X_val_i = X_val
    L = X_val_i.shape[0]
    # iterate over bond types 
    for b in bond_types:
        # unpacking hyperparameters and pseudo inputs
        M_bond = bond_dict[b]['M']
        X_pseudo_bond = bond_dict[b]['X_pseudo']
        delta_bond = bond_dict[b]['delta']
        theta_bond = bond_dict[b]['theta']
        beta_bond = bond_dict[b]['beta']
        cutoff_bond = bond_dict[b]['cutoff']
        d_bond = bond_dict[b]['d']
        # iterate over pseudo inputs of current bond type
        for j in range(M_bond):
            X_pseudo_bond_j = X_pseudo_bond[j]
            temp = 0.
            for l in range(L):
                # skip nonequal bond-types (they have zero covariance)
                if B_val_i[l] != b:
                    continue
                temp += k.evaluate(X_val_i[l], X_pseudo_bond_j,
                                                delta_bond,
                                                theta_bond,
                                                beta_bond,
                                                cutoff_bond,
                                                d_bond)
            k_star[j+offset] = temp

        offset += M_bond

    return np.dot(k_star, alpha)

def posterior_variance(Kernel k,
                       bond,
                       bond_types,
                       bond_dict,
                       int M,
                       np.ndarray[np.float64_t, ndim=2] inv_K_M,
                       np.ndarray[np.float64_t, ndim=2] inv_Q_M,
                       np.ndarray[np.float64_t, ndim=1] X_star):

    cdef np.ndarray[np.float64_t, ndim=1] k_star = np.zeros(M, dtype=np.float64)
    cdef int P = X_star.shape[0]
    cdef np.ndarray[np.float64_t, ndim=1] result = np.zeros(P, dtype=np.float64)
    cdef np.ndarray[np.float64_t, ndim=1] X_pseudo_bond
    cdef Py_ssize_t i, j
    
    cdef np.ndarray[np.float64_t, ndim=2] A = inv_K_M - inv_Q_M
    

    # determine start index for the M_bond non-zero elements of k_star
    index_of_bond = bond_types.index(bond)
    offset = 0
    for i in range(index_of_bond):
        offset += bond_dict[bond_types[i]]['M']

    # unpacking hyperparameters and pseudo inputs
    M_bond = bond_dict[bond]['M']
    X_pseudo_bond = bond_dict[bond]['X_pseudo']
    delta_bond = bond_dict[bond]['delta']
    theta_bond = bond_dict[bond]['theta']
    beta_bond = bond_dict[bond]['beta']
    cutoff_bond = bond_dict[bond]['cutoff']
    d_bond = bond_dict[bond]['d']

    # calculate result
    for i in range(P):
        xi = X_star[i]
        for j in range(M_bond):
            k_star[j+offset] = k.evaluate(xi, X_pseudo_bond[j],
                                          delta_bond, 
                                          theta_bond, 
                                          beta_bond, 
                                          cutoff_bond, 
                                          d_bond)
        result[i] = (k.evaluate(xi, xi,
                                delta_bond,
                                theta_bond,
                                beta_bond,
                                cutoff_bond,
                                d_bond) - np.dot(k_star, np.dot(A, k_star))) 
    
    return result
    
